knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(kableExtra)
library(mgcv)
library(ggpmisc)
library(mgcViz)
library(viridis)

source('scripts/TSS_mclapply_ie_stn_month.R')

load.Rdata <- function( filename, objname )
{
  d1 <- load( filename )
  eval( parse( text=paste( objname,  "<<- ", d1  ) ) )
}

dat<- as.data.frame(read.delim("data/cal_vert_data_std_light_bloom_for_GAMs.txt"))

month_assign_cfin=data.frame(Month=1:12, TAXON="Cfin") %>% mutate(season=ifelse(Month >6| Month <3, "D", "A"))
month_assign_chyp=data.frame(Month=1:12, TAXON="Chyp") %>% mutate(season=ifelse(Month >3, "D", "A"))
month_assign_cglac=data.frame(Month=1:12, TAXON="Cglac") %>% mutate(season=ifelse(Month %in% c(4,5), "A", "D"))

month_assign<- rbind(month_assign_cfin,month_assign_chyp,month_assign_cglac)



dat<- left_join(dat %>%  mutate(fMonth=as.factor(Month),
                                             cum_n_m2 = ind.m2_total * pcum), month_assign) # for simulation



theme_cl<- theme_few()+theme(axis.title.x = element_text(size=14, colour="black"),
                             axis.text.x = element_text(size=13, colour='black'),
                             axis.title.y = element_text(size=14, colour="black"),
                             axis.text.y = element_text(size=13, colour='black'),
                             legend.title=element_blank(),
                             panel.border = element_rect(colour = "black"),
                             axis.ticks = element_line(colour="black"),
                             strip.text.x = element_text(size=14, face="bold"),
                             strip.text.y = element_text(size=14, face="bold"))

reducing eps in gam formula with betar distribution has a positive effect on the normality of random effects.

What we need to know about the data from the EcoMon survey

criteria for correction: 1) % of depth sampled < 95% (tow_depth/sta_depth) 2) sta_depth-tow_depth > 10. if sta_depth-tow_depth < 10, then considered as sampled properly.

#to be replaced by Ben package?
library(tidyxl)
ecomon<- readxl::read_excel(sheet="Data","C:/LEHOUX/Projects/Données_copépodes/ACCASP2020/Calfin_Thru_12_30_2019_10m2.xlsx")%>%  mutate(perc_sampled= tow_depth/sta_depth *100,
                                                       to_correct=ifelse(perc_sampled < 95 & (sta_depth-tow_depth) >10, "Y", "N"),
                   Zcat = as.factor(ifelse(sta_depth < 200, "0-200",
                                           ifelse(sta_depth >=200 & sta_depth < 300, "200-300",
                                                  ifelse(sta_depth >=300 & sta_depth < 500, "300-500",
                                                         ifelse(sta_depth>=500 & sta_depth < 1000, "500-1000", 
                                                             ifelse(sta_depth>=1000, "1000+" , NA)))))))


ecomon %>%  group_by(to_correct,Zcat) %>%  tally()
## # A tibble: 8 × 3
## # Groups:   to_correct [2]
##   to_correct Zcat         n
##   <chr>      <fct>    <int>
## 1 N          0-200    20502
## 2 N          200-300    627
## 3 N          300-500      3
## 4 Y          0-200     1403
## 5 Y          1000+        5
## 6 Y          200-300   1563
## 7 Y          300-500    269
## 8 Y          500-1000    15

Cfin CVI

load("gamms/Cfin_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      73.65314  -980.0647
## mod4      47.62006 -2209.9254
## modseason 60.45941  -950.3503
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F)

Profiles

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod3),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod3, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Month as continuous

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod4),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Season

dattax_season <-  dattax_season %>%  mutate(residuals=residuals(modseason),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

dattax_season$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Vis.gam

Month as factor

print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")

print(plot(getViz(mod3), all.terms=T), pages=1)

new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
                          Zstation = seq(50,500,10),
                          percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
                     as.data.frame(predict(mod3, newdata = new_data,
                                           se.fit = TRUE)))

ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
                    colour = fMonth)) +
    geom_line() + theme_cl+
    facet_wrap(~ fMonth) 

#re<- as.data.frame(predict(mod3, type="terms")) %>%  select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+ 
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

datsp <-  dat %>%  filter(TAXON_STADE=="Cfin_CVI")

Month as factor

dat200 <-  datsp %>%  dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>%   dplyr::slice(which.max(ZMax)) 

dat_zstation <-  datsp %>% dplyr::group_by(ID) %>%   slice(which.max(ZMax))

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Month as continuous

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Season

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Cfin CV

load("gamms/Cfin_CVgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CVgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      237.2257  -4282.196
## mod4      149.0354 -11564.548
## modseason 157.3500  -4231.050
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F)

Profiles

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod3),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod3, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Month as continuous

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod4),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Season

dattax_season <-  dattax_season %>%  mutate(residuals=residuals(modseason),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

dattax_season$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Vis.gam

Month as factor

print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")

print(plot(getViz(mod3), all.terms=T), pages=1)

new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
                          Zstation = seq(50,500,10),
                          percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
                     as.data.frame(predict(mod3, newdata = new_data,
                                           se.fit = TRUE)))

ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
                    colour = fMonth)) +
    geom_line() + theme_cl+
    facet_wrap(~ fMonth) 

#re<- as.data.frame(predict(mod3, type="terms")) %>%  select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+ 
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

datsp <-  dat %>%  filter(TAXON_STADE=="Cfin_CV")

Month as factor

dat200 <-  datsp %>%  dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>%   dplyr::slice(which.max(ZMax)) 

dat_zstation <-  datsp %>% dplyr::group_by(ID) %>%   slice(which.max(ZMax))

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Month as continuous

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Season

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Cfin CIV

load("gamms/Cfin_CIVgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CIVgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      185.3930  -3763.047
## mod4      118.6907 -10792.087
## modseason 126.7773  -3780.707
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F)

Profiles

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod3),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod3, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Month as continuous

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod4),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Season

dattax_season <-  dattax_season %>%  mutate(residuals=residuals(modseason),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

dattax_season$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Vis.gam

Month as factor

print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")

print(plot(getViz(mod3), all.terms=T), pages=1)

new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
                          Zstation = seq(50,500,10),
                          percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
                     as.data.frame(predict(mod3, newdata = new_data,
                                           se.fit = TRUE)))

ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
                    colour = fMonth)) +
    geom_line() + theme_cl+
    facet_wrap(~ fMonth) 

#re<- as.data.frame(predict(mod3, type="terms")) %>%  select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+ 
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

datsp <-  dat %>%  filter(TAXON_STADE=="Cfin_CIV")

Month as factor

dat200 <-  datsp %>%  dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>%   dplyr::slice(which.max(ZMax)) 

dat_zstation <-  datsp %>% dplyr::group_by(ID) %>%   slice(which.max(ZMax))

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Month as continuous

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Season

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Cfin CIV-CVI

load("gamms/Cfin_CIV_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CIV_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      276.4344  -4488.609
## mod4      153.7038 -10592.331
## modseason 194.2298  -4322.528
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F) 
tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

  kbl(tss_table) %>%
  kable_paper(full_width = F)
ObsPrev_abund Tresh_opt_e_i TSS_e_i TSSsd_e_i AUC_e_i corr_r_e_i p.corr_e_i intercept_e_i slope_e_i r2_e_i r2sd_e_i
0.65 NA/NA NA/NA NA/NA NA/NA 0.89/0.89 0/0 0.12/0.16 0.79/0.73 0.77/0.76 0.01/0
0.65 NA/NA NA/NA NA/NA NA/NA 0.89/0.86 0/0 0.16/0.2 0.73/0.67 0.74/0.7 0.05/0
0.65 NA/NA NA/NA NA/NA NA/NA 0.85/0.86 0/0 0.23/0.21 0.66/0.64 0.67/0.69 0.02/0

Profiles

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod3),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod3, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Month as continuous

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod4),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Season

dattax_season <-  dattax_season %>%  mutate(residuals=residuals(modseason),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

dattax_season$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Vis.gam

Month as factor

print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")

print(plot(getViz(mod3), all.terms=T), pages=1)

new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
                          Zstation = seq(50,500,10),
                          percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
                     as.data.frame(predict(mod3, newdata = new_data,
                                           se.fit = TRUE)))

ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
                    colour = fMonth)) +
    geom_line() + theme_cl+
    facet_wrap(~ fMonth) 

#re<- as.data.frame(predict(mod3, type="terms")) %>%  select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+ 
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

datsp <-  dat %>%  filter(TAXON_STADE=="Cfin_CIV_CVI")

Month as factor

dat200 <-  datsp %>%  dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>%   dplyr::slice(which.max(ZMax)) 

dat_zstation <-  datsp %>% dplyr::group_by(ID) %>%   slice(which.max(ZMax))

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Month as continuous

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Season

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Cfin CV-CVI

load("gamms/Cfin_CV_CVIgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cfin_CV_CVIgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      278.2126  -4480.242
## mod4      153.4178 -10597.444
## modseason 194.1278  -4329.660
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F) 
tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

  kbl(tss_table) %>%
  kable_paper(full_width = F)
ObsPrev_abund Tresh_opt_e_i TSS_e_i TSSsd_e_i AUC_e_i corr_r_e_i p.corr_e_i intercept_e_i slope_e_i r2_e_i r2sd_e_i
0.65 NA/NA NA/NA NA/NA NA/NA 0.89/0.89 0/0 0.12/0.16 0.79/0.73 0.76/0.76 0.01/0
0.65 NA/NA NA/NA NA/NA NA/NA 0.89/0.86 0/0 0.16/0.2 0.73/0.67 0.74/0.7 0.04/0
0.65 NA/NA NA/NA NA/NA NA/NA 0.85/0.86 0/0 0.22/0.21 0.66/0.64 0.68/0.7 0.02/0

Profiles

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod3),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod3, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod3, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Month as continuous

dattax_month <-  dattax_month %>%  mutate(residuals=residuals(mod4),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_month, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point()+
  geom_smooth()

 dattax_month$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_month, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_month, aes(y=predict(mod4, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Season

dattax_season <-  dattax_season %>%  mutate(residuals=residuals(modseason),
                                          cat_perc= ifelse(percZ_stn < 0.25, "<25",
                                                           ifelse(percZ_stn >=0.25 & percZ_stn < 0.50, "25-50",
                                                                  ifelse(percZ_stn >= 0.50 & percZ_stn < 0.75, "50-75", ">75"))))


ggplot(dattax_season, aes(y=residuals, x=Month))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=LONGITUDE))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Year))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

ggplot(dattax_season, aes(y=residuals, x=Zstation))+ facet_wrap(~cat_perc)+
  geom_point(aes(col=season))+
  geom_smooth()

dattax_season$linpred<- as.vector(predict(mod4, type="link"))

ggplot(dattax_season, aes(y=residuals, x=linpred))+ ggtitle("Homogeneity of residuals")+
  geom_point()+
  geom_smooth()

ggplot(dattax_season, aes(y=predict(modseason, type="response"), x=pcum))+ ggtitle("Prediction accuracy")+
  geom_point()+
  geom_smooth()

Vis.gam

Month as factor

print(plot(getViz(mod3), all.terms=T, select=ncol(predict(mod3, type="terms"))))

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="3"), main="Month3")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="8"), main="Month8")

vis.gam(mod3, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(fMonth="10"), main="Month10")

print(plot(getViz(mod3), all.terms=T), pages=1)

new_data <- tidyr::expand(dattax_month, nesting(ID,fMonth),
                          Zstation = seq(50,500,10),
                          percZ_stn=0.5)
m1_pred <- bind_cols(new_data,
                     as.data.frame(predict(mod3, newdata = new_data,
                                           se.fit = TRUE)))

ggplot(m1_pred, aes(x = Zstation, y = fit, group = ID,
                    colour = fMonth)) +
    geom_line() + theme_cl+
    facet_wrap(~ fMonth) 

#re<- as.data.frame(predict(mod3, type="terms")) %>%  select(`s(ID)`)
#dattax_month$re<- re$`s(ID)`
#dattax_month$residuals<-residuals(mod3)#

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=ind.m2_total))+ scale_color_viridis(option="turbo",trans="log10") +theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=LONGITUDE)) +
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Zstation))+ scale_color_viridis(option="turbo",trans="log10")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=Year))+ 
#scale_color_viridis(option="turbo")+theme_cl

#ggplot(dattax_month, aes(x=residuals, y=re))+geom_point(aes(col=fMonth)) +
#scale_color_viridis_d(option="turbo")+theme_cl

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

datsp <-  dat %>%  filter(TAXON_STADE=="Cfin_CV_CVI")

Month as factor

dat200 <-  datsp %>%  dplyr::filter(ZMax < 200 & Zstation > 200) %>% dplyr::group_by(ID) %>%   dplyr::slice(which.max(ZMax)) 

dat_zstation <-  datsp %>% dplyr::group_by(ID) %>%   slice(which.max(ZMax))

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod3, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Month as continuous

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(mod4, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Season

dat200$corrected_bioness =dat200$cum_n_m2/as.vector(predict(modseason, dat200, type="response", exclude ='s(ID)', newdata.guaranteed = T))
dat_corr<- left_join(dat_zstation, dat200[,c("ID", "corrected_bioness")])

dat_corr <-  dat_corr %>%  mutate(Zcat = as.factor(ifelse(Zstation < 300, "< 300",
                                                  ifelse(Zstation >=300 & Zstation < 500, "300-500",
                                                         ifelse(Zstation>=500 & Zstation < 1000, "500-1000", 
                                                             ifelse(Zstation>=1000, ">1000" , NA))))) , Zcat=recode_factor(Zcat, `<300` ="<300", `300-500`="300-500",`500-1000`="500-1000",`>1000`=">1000"))

ggplot(data=dat_corr, aes(y=corrected_bioness, x=cum_n_m2))+
  geom_point(aes(col=Zcat))+ #montre pas de lien avec factor year
  geom_smooth(method="lm", se=F)+
  geom_abline(slope=1, intercept=0)+
  theme_few()+
  stat_poly_eq(formula = y~x,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +   
  scale_y_continuous(trans="log10", name="Simulated corrected integrated abundance")+
  scale_x_continuous(trans="log10", name="Observed intergrated abundance")

Chyp CIV-CVI

load("gamms/Chyp_CIV_VIgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Chyp_CIV_VIgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Chyp_CIV_VIgammsall_remove01.RData")

Table

AIC(mod3,mod4, modseason)
##                 df       AIC
## mod3      52.83426 -1135.271
## mod4      25.05714 -3068.091
## modseason 29.95230 -1126.619
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F) 
AIC(mod3,mod4, modseason, modall)
##                 df       AIC
## mod3      52.83426 -1135.271
## mod4      25.05714 -3068.091
## modseason 29.95230 -1126.619
## modall    26.31336 -1131.882
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

 # kbl(tss_table) %>%
  #kable_paper(full_width = F)

Profiles

No interaction

no int

Month as factor

failed at the prediction level

Month as continuous

month as continuous

Season

Season

Residuals patterns

Month as factor

Month as continuous

Season

Vis.gam

No interactions

print(plot(getViz(modall), all.terms=T, select=ncol(predict(modall, type="terms"))))

vis.gam(modall, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray")

print(plot(getViz(modall), all.terms=T), pages=1)

Month as factor

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

continuous

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

Month as factor

Month as continuous

Season

Cgla CIV-CVI

load("gamms/Cglac_CIV_VIgamms3_4_remove01.RData"); dattax_month <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cglac_CIV_VIgammsseason_remove01.RData"); dattax_season <- dattax %>%  mutate(fMonth=as.factor(Month))
load("gamms/Cglac_CIV_VIgammsall_remove01.RData")

Table

AIC(mod3,mod4, modseason)
##                 df        AIC
## mod3      25.97025  -617.0373
## mod4      16.88213 -1842.7810
## modseason 21.41561  -598.0254
#doesn't work don't know why?? all factor levels are there. doesn't even wirk with predict dataC 
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.0000001)", nrem=4,Sp="Cfin_CVI")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

#  kbl(tss_table) %>%
#  kable_paper(full_width = F) 
AIC(mod3,mod4, modseason, modall)
##                 df        AIC
## mod3      25.97025  -617.0373
## mod4      16.88213 -1842.7810
## modseason 21.41561  -598.0254
## modall    20.63003  -597.3538
#tss_tablemod3<- TSS_mclapply_ie_stn(mcomp=mod3, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemod4<- TSS_mclapply_ie_stn(mcomp=mod4, Data=dattax_month, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_tablemodseason<- TSS_mclapply_ie_stn(mcomp=modseason, Data=dattax_season, gam_type="standard", fam="betar(eps=0.001)", nrem=4,Sp="pcum")  
#tss_table<- rbind(tss_tablemod3,tss_tablemod4,tss_tablemodseason)

 # kbl(tss_table) %>%
  #kable_paper(full_width = F)

Profiles

No interaction

month as factor

Month as factor

month as factor ### Month as continuous month as continuous ### Season month as season

Residuals patterns

Month as factor

Month as continuous

Season

Vis.gam

No interactions

print(plot(getViz(modall), all.terms=T, select=ncol(predict(modall, type="terms"))))

vis.gam(modall, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray")

print(plot(getViz(modall), all.terms=T), pages=1)

Month as factor

Month as continuous

print(plot(getViz(mod4), all.terms=T, select=ncol(predict(mod4, type="terms"))))

season

Season

print(plot(getViz(modseason), all.terms=T, select=ncol(predict(modseason, type="terms"))))

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="D"), main="Diapause")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

vis.gam(modseason, view=c("Zstation", "percZ_stn"),plot.type="contour", type="response", color="gray",cond=list(season="A"), main="Active")

print(plot(getViz(modseason), all.terms=T), pages=1)

Simulation of bioness data

Month as factor

Month as continuous

Season